Method for analyzing multi phase and heat flow of fluids in reservoir and recording media therefor

ABSTRACT

A method for analyzing multi phase and heat flow in a reservoir is disclosed. 
     The method comprises a stage for specifying fracture surfaces in a reservoir as multiple three-dimensional coordinates, and specifying a three-dimensional coordinates of a intersecting line which intersects between the fracture surfaces to grasp connectivity between the fracture surfaces, a stage for generating a mathematical equation having main variables of pressure, saturation and temperature in the cell from at least two or more of fluids among oil, water, gas that are moving within the fracture network and each flow model of fluids with heat, and a stage for grasping the flow of fluids and heat within the fracture network in the reservoir by activating the flow model by specifying initial values of pressure, saturation and temperature in each cell, and specifying production condition or injection condition of a drilling well connected to the fracture network.

TECHNICAL FIELD

The present invention relates to both resources development technologies for development of resources such as oil, gas endowed in a reservoir and to environmental technologies such as CCS for storage of carbon dioxide in the ground, and more particularly to a method of grasping and estimating a flow of a multi phase fluid including heat in an underground reservoir in accordance with variable conditions.

BACKGROUND ART

The importance of developing unconventional gas such as shale gas, tight gas and coal bed methane has been increasing among oil development companies due to price rise of conventional resources such as oil and gas along with sharply increasing energy demands over a past decade. For example, shale gas is being actively produced in North American regions such as Barnett, Haynesville, Woodford and Engleford. It is foreseen that shale gas will replace a great portion in fossil fuel in the future, and therefore resource development companies, investment companies, and research companies show big interests in the shale gas.

A reason why a commercial production of the shale gas has become possible can be found in technical advancement which enabled a horizontal drilling which is an essential means in shale gas development being carried out in an economical way and in the market situation of increasing energy demands or oil price rise. The horizontal drilling is a technology for digging through a reservoir in a horizontal direction along the sweet spots in the reservoir.

FIG. 1 is a schematic view representing endowment characteristics and drilling development of conventional gas and unconventional gas such as shale gas, and FIG. 2 is a perspective schematic view of a shale gas field as viewed from above.

Referring to FIGS. 1 and 2, the shale gas is characterized to be endowed in a tight shale reservoir having very small porosity and permeability while the conventional resources such as oil and gas are endowed in reservoir rocks having big porosity and permeability. Further, the shale gas is characterized to be distributed in a wide horizontal region extended along the shale while the conventional resources such as oil or gas are intensively endowed in a certain region.

Hence, a horizontal drilling in multiple directions from one vertical production well is to be carried out as shown in FIG. 2, and fractures are to be formed by consecutively fracturing the reservoir rock at a short interval along the horizontal drilling section. The shale gas is introduced into a horizontal drilling well along the fractures and then produced onto the ground through a vertical drilling well.

As described above, the shale gas is characterized to move through the interconnected fractures (by hydraulic fracturing or natural formed) rather than flowing through pores since the reservoir rocks have very low porosity in the case of shale gas while the oil flows through pores of the reservoir rocks when the pressure inside the reservoir rocks is open by drilling a vertical drilling well since the reservoir rocks have very high porosity and permeability in the case of the conventional resources such as oil or natural gas.

In developing oil, natural gas as well as shale gas, economic validity is firstly studied on the basis of data retrieved from geophysical survey since it is an important matter to predict how the fluid will behave in a reservoir under a specific condition applied when a drilling well is formed. However, most of fluid flow models conventionally developed are aimed at conventional resources like oil or natural gas. That is, the existing analysis methods are designed to simulate movement of fluids on the basic premise of fluids moving through pores in the reservoir, and therefore there has been a limitation that existing analysis methods are not adapted to the flow model for fluids, such as shale gas, moving through the fractures.

Further, there is another limitation that the part of heat flow is missed from some programs that have developed for modeling the flow of unconventional resources such as shale gas. In the case of conventional oil or natural gas, the fluid flow represents far different aspect by temperature condition associated with the heat flow. Therefore, it is presumed that fluids such as shale gas moving through fractures also represent different flow characteristic in the case of the heat flow existence. However, many limitations have been exposed in reflecting real phenomenon with the programs modeled the flow of shale gas due to lack of heat flow consideration.

DISCLOSURE Technical Problem

The present invention has been devised to resolve limitations described above, and is to provide a method for analyzing multi phase and heat flow of fluids in a reservoir in which the multi phase flow occurs through an interconnected fracture network and the heat flow is particularly reflected therein in order to model the real flow of unconventional gas like shale gas to be very similar to the real flow.

Advantageous Effects

In the present invention, a discrete fracture network model adapted to model fluids that are moving through fractures is used, and oil, water, gas and heat are modeled as fluids moving along the fractures in the reservoir in order to copy the multi phase flow. The significance of the present invention is relevant for provision of a base enabling the prediction of various real aspects of the shale gas flow in the reservoir by such a method.

Above all, the present invention enables a more reliable modeling for the fluid flow in the shale gas area by adding the heat flow to an existing discrete fracture network flow model since there has not been enough research carried out for the part of heat and temperature in the case of multi phase flow in the existing discrete fracture network flow model despite the temperature condition acting as a very important factor in the multi phase flow.

DESCRIPTION OF DRAWINGS

FIG. 1 is a schematic diagram showing endowment characteristics and drilling development of conventional gas and unconventional gas such as shale gas.

FIG. 2 is a schematic diagram illustrating a vertical drilling well and a horizontal drilling well of a shale gas pad.

FIG. 3 is a diagrammatic flow chart illustrating a method for analyzing a multi phase and a heat flow of fluids in a reservoir in accordance with an embodiment of the present invention.

FIG. 4 is a perspective view illustrating a reservoir partitioned into cells in a dual porosity reservoir model.

FIG. 5 is a perspective view illustrating a fracture surface and an intersecting reference line characterized in an embodiment of the present invention.

FIG. 6 is a perspective view illustrating a cell on a fracture surface partitioned along an intersecting line between fracture surfaces.

BEST MODE Mode for Invention

A more detailed description for a method analyzing multi phase and heat flow of fluids in a reservoir in accordance with an embodiment of the present invention will be described below. The present invention relates to a time series method, and is implemented by a computer-readable program. Accordingly, the present invention is executed by activating a program installed in a computer or recorded in a recording media through a computer.

FIG. 3 is diagrammatic flow chart of a method (referred to as “multi phase flow analysis method” hereinafter) for analyzing multi phase and heat flow of fluids in a reservoir in accordance with an embodiment of the present invention.

Referring to FIG. 3 the multi phase flow analysis method begins with selection of a site (e.g. a shale gas reservoir) for analyzing fluids flow, and then grasping fractures in the reservoir that is a three dimensional space.

In general, conventional resources such as oil or natural gas move through pores within the reservoir. Accordingly, the entire reservoir is partitioned into multiple cubic-shaped cells by a dual porosity reservoir model, and the fluid flow is copied by applying conditions such as porosity, permeability and fluids pressure. However, the fluids in a reservoir such as a shale gas reservoir having a tight lithology are characterized by moving through fractures rather than pores due to very low porosity and permeability of the reservoir. Accordingly, it is preferred to implement a discrete fracture network flow model instead of implementing a dual porosity reservoir model in the case of modeling reservoirs composed of shale.

In order to implement the discrete fracture network model, the fractures within the reservoir are firstly figured out. The targeted fractures in the case of shale gas are the natural fractures in the reservoir as well as the continuous fractures within a horizontal drilling section formed by a hydraulic fracturing. In the present invention, the fracture is assumed as a surface in three-dimensional space. Accordingly, the fracture surfaces (S1 and S2) can be possibly specified by four sets of three-dimensional coordinates (S1: A, B, C, D, S2: O, P, Q, R) in the case of a fracture on a rectangular surface as shown for example in FIG. 5. The fracture surfaces are, however, not the surfaces having a mathematical meaning, and the fracture surfaces have thickness and volume as well. As thickness is constant, the fracture can be specified as a shape having thickness on sets of coordinates as described above.

As described above, an intersecting line between respective fracture surfaces is grasped after multiple fracture surfaces in reservoir are specified. The intersecting line between fracture surfaces is very important in fluid flow. When a drilling well is formed in a reservoir, the pressure equilibrium is broken so that fluids move to the drilling well along the interconnected fracture surfaces from the original location. However, a fluid on an independent fracture surface, among fracture surfaces, which is unconnected with other fracture surfaces, can not continuously move since the fluid is limited in movement to be within the independent fracture surface.

Therefore, a possible production amount of the shale gas can only be figured out in the case that the shale gas is led from the original location to production well finally by moving along the fracture surfaces, i.e. the fracture surfaces are continuously led from the beginning point to the production well. That is, the connectivity of the fractures is a very important factor in terms of the production amount of oil or gas.

As the fracture surfaces are specified as coordinates in the three-dimensional space, the intersecting line (L) between the fracture surfaces (S1, S2) can be easily and mathematically found out as illustrated in FIG. 4. And the intersecting line is again specified as two coordinates in three-dimensional space.

As described above, when respective intersecting lines produced by multiple fracture surfaces intersecting each other are all figured out, each of the fracture surfaces is partitioned into multiple cells using Finite Different Method (FDM), and a fractures network is finally formed. Cells on each of the fracture surfaces may be formed to be identical in size and shape or be combined to have various size and shape. What is important is that, as shown in FIG. 6, the cells (1, 2, 3 . . . ) are partitioned by the intersecting line (L) between the fracture surfaces. That is, the intersecting line (L) does not pass through the inside of cells. Conversely, the partitioning of the cells is performed by a grid partitioning with the intersecting line (L) as a boundary line.

As described above, the fracture network in the reservoir is produced when the fracture surfaces and the intersecting line in the reservoir are specified, and cells are formed by partitioning each of the fracture surfaces. And, conditions, i.e., fluids pressure and fluids saturation rate to determine fluids flow are introduced to each cells. More particularly, this embodiment relates to copy multi phase flow of oil, gas, water and heat flowing together along the fractures and therefore pressure of the fluids (or partial pressure), relative saturation rate of each fluid, and temperature are predetermined in each cells. And, when the fluid flows, these values continuously change with the lapse of time and copy the flow of the fluid. For example, oil flow direction is determined by oil pressure (high or low) in two adjacent cells, and heat transfer direction is determined by temperature difference. And, as for oil production, the production amount can be figured out by tracing the amount of oil movement with the lapse of time from a cell (higher in pressure than a drilling well) adjacent to a drilling well cell to the drilling well cell to which the drilling well is connected.

That is, the condition of the cell with the lapse of time, i.e., pressure of each fluid, saturation and temperature change with the lapse of time, and these changes subsequently allow fluid and heat to flow, and such a change that is mathematically modeled is the flow model.

First of all, the flow model of oil is represented by non-linear partial differential equations as shown in Equation 1 below, and the flow model of water is shown in Equation 2.

$\begin{matrix} {{{{\frac{\partial}{\partial L}\left\lbrack {T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; W} + Q_{o} + Q_{Io}} = {\left( {\varphi \; V_{f}} \right)\frac{\partial}{\partial t}\left( \frac{S_{o}}{B_{o}} \right)}} & {{Equation}\mspace{14mu} 1} \\ {{{{\frac{\partial}{\partial L}\left\lbrack {T_{wL}\frac{\partial}{\partial L}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {T_{wW}\frac{\partial}{\partial W}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; W} + Q_{w} + Q_{Iw}} = {\left( {\varphi \; V_{f}} \right)\frac{\partial}{\partial t}\left( \frac{S_{w}}{B_{w}} \right)}} & {{Equation}\mspace{14mu} 2} \end{matrix}$

where L is a cell length, W is a cell width, T_(L) is a longitudinal conductivity of fluid in cell, T_(W) is width directional conductivity of fluid in cell, p is a fluid pressure in cell, γ is a cell density, z is a cell depth, Q is an amount of fluid produced in or injected from the drilling well, Q_(I) is a flow amount of fluid on intersecting line between fracture surfaces, φ is a cell porosity, V_(f) is a volume of a fracture surface, t is time, S is a saturation rate of fluid in cell, B is a volume coefficient (volume) of fluid in cell, subscript o represents oil, and subscript w represents water.

The above Equation 1 and 2 are identical except for kinds of fluids. That is, there are four articles on the left side and one article on the right side.

The first two of the four articles on the left side represent flow amount of the fluid in each of longitudinal direction and width direction. And Q of the latter two of four articles on the left side is production amount or injection amount at the drilling well which can be defined as flow amount of the fluid within the cell, and Q_(I) is flow amount of the fluid on the intersecting line between fracture surfaces. But, as a value of Q exists in a cell through which a drilling well passes, and therefore the value of Q becomes zero in most of other cells, and likewise, Q_(I) is flow amount of the fluid on intersecting line between fracture surfaces, the value of Q_(I) exists only in cells in contact with the intersecting line. Q and Q_(I) are for a case of special flow of the fluid existing in a cell under special conditions (either a drilling well passes through or the cell is in contact with intersecting line). That is, the left sides of both Equation 1 and 2 are defined as the normal flow amount from longitudinal direction and width direction added with flow amount of the fluids under special conditions such as being passed through by the drilling well or the intersecting line.

And, the right side of Equation 1 and 2 represent the variation of saturation (S) inside the cell along the lapse of time (T). For reference, the porosity (φ) on the right side is assumed as constant such as volume of the fracture (Vf) since the porosity is on the premise of a case that fillers are stacked on the fracture surfaces.

Finally, the flow model of Equation 1 and 2 can represent the relation between a variation of flow amount and a variation of saturation within the cell, and a main variable on the left side is the pressure (p) and a main variable on the right side is the saturation (S). That is, the flow model is defined by the variation of the pressure and the saturation. For reference, conductivity (T) or formation volume factor (B), Q, and Q_(I) are dependent variables being set by the main variables of the pressure and the saturation that are separately defined.

Meanwhile, like the flow model of oil and water, a gas flow model is a non-linear partial differential equation as described in Equation 3 below.

$\begin{matrix} {{{\frac{\partial}{\partial L}\left\lbrack {T_{gL}\frac{\partial}{\partial L}\left( {p_{g} - {\gamma_{g}z}} \right)} \right\rbrack} + {R_{so}T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {{T_{gW}\frac{\partial}{\partial W}\left( {p_{g} - {\gamma_{g}z}} \right)} + {R_{so}T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)}} \right\rbrack}\Delta \; W} + Q_{g} + {R_{so}Q_{o}} + Q_{Ig} + {R_{so}Q_{Io}}} = {\left( {\varphi \; V_{f}} \right)\frac{\partial}{\partial t}\left( {\frac{S_{g}}{B_{g}} + \frac{R_{so}S_{o}}{B_{o}}} \right)}} & {{Equation}\mspace{14mu} 3} \end{matrix}$

where symbols herein are identical to those of previous Equation 1 and 2 except for R_(s) which represents a gas dissolution ratio in other fluid, and subscript g represents the gas.

Examining the gas flow model, the gas flow model is identical to those of oil or water flow model in a sense that articles for variation of longitudinal direction flow amount and width direction flow amount of a cell, injection amount at a drilling well and flow amount at an intersecting line are included. However, there is an additional article to be added in comparison to oil or water flow models. That is, the amount of gas dissolved in oil has to be further included. Accordingly, in gas flow model, the article of oil flow amount multiplied by dissolution rate (Rs) is further added together with gas flow amount itself. Likewise, the right side is further added with an article of oil saturation rate multiplied by dissolution rate together with a relative saturation rate of gas itself. The rest articles will not be described for these are identical to those of oil and water flow models.

Meanwhile, there have been very few cases that the heat flow is considered in a discrete fracture network model in conventional arts, however, the fluidity of fluids is represented in a wide variety of aspects depending on the temperature (heat flow) condition in a real shale gas field, and therefore the heat flow model needs to be included in the reservoir fracturing model. Hence, the present invention has reached to inclusion of the heat flow model in the reservoir fracturing model as shown in Equation 4 below.

$\begin{matrix} {{{{\frac{\partial}{\partial L}\left\lbrack {H_{o}T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {H_{w}T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; W} + {{\frac{\partial}{\partial L}\left\lbrack {{H_{g}T_{gL}\frac{\partial}{\partial L}\left( {p_{g} - {\gamma_{g}z}} \right)} + {H_{g}T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)}} \right\rbrack}\Delta \; {L++}{\frac{\partial}{\partial W}\left\lbrack {{H_{g}T_{gW}\frac{\partial}{\partial W}\left( {p_{g} - {\gamma_{g}z}} \right)} + {H_{o}R_{so}T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)}} \right\rbrack}\Delta \; W} + {{\frac{\partial}{\partial L}\left\lbrack {H_{w}T_{wL}\frac{\partial}{\partial L}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {H_{w}T_{wW}\frac{\partial}{\partial W}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; W} + {H_{o}\left( {Q_{o} + Q_{Io}} \right)} + {H_{g}\left\lbrack {Q_{g} + Q_{Ig} + {R_{so}\left( {Q_{o} + Q_{Io}} \right)}} \right\rbrack} + {H_{w}\left( {Q_{w} + Q_{Iw}} \right)} + {\frac{\partial}{\partial L}\left( {T_{h}\frac{\partial\theta}{\partial L}} \right)} + {\frac{\partial}{\partial W}\left( {T_{h}\frac{\partial\theta}{\partial W}} \right)}} = {\varphi \; V_{f}\frac{\partial}{\partial t}(U)}} & {{Equation}\mspace{14mu} 4} \end{matrix}$

where symbols herein are identical to those of previous Equation 1, 2, and 3 while an additional H represents an enthalpy of a fluid, U represents an inner energy in cell, θ represents a temperature, and subscript h represents heat.

The heat flow model has an article for the temperature variation amount (as appeared on the last line of Equation 4) along the longitudinal direction and the width direction like other fluids previously described.

Further, the heat movement amount transferred during the fluids of oil, water, gas move is also included. That is, the first line of Equation 4 is for the oil flow amount multiplied by the enthalpy (H), the second and third lines are for the gas flow amount multiplied by enthalpy, and the fourth line is for water flow amount multiplied by enthalpy. And, the fifth line is for the sum of movement of oil, water, and gas multiplied by enthalpy at the drilling well and on the intersecting line. That is, the first through the fifth lines are for flow amount of each of oil, gas, and water multiplied by enthalpy, the heat movement amount can be defined as the fluids flow amount multiplied by the enthalpy since the enthalpy is the heat content (or heat transfer amount) of the fluids.

That is, the heat flow model is defined to include both the temperature variation amount along the longitudinal direction and the width direction and the heat amount transferred by the fluid. And the right side of the Equation 4 is defined as the variation of internal energy (U) within the cell along the lapse of time. The internal energy is a function of the temperature in which the temperature variation in the cell by the heat flow is reflected.

For reference, the Equations 1-4 are all non-linear partial differential equations that can be solved by finding an approximate value. That is, when the value for the partial differential equation is found, the value is verified whether or not to be in an error range, and then a process of finding the approximate value repeats by inputting the value found back as an initial value if the value found is verified to be out of the error range, and finally the approximate value in the error range is found. Mathematically, the partial differential equation is transformed to an ordinary differential equation and the value can be found.

In the present invention, the values of each fluid flow model are pressure, saturation rate and the temperature. The result of the flow model may be finally represented as variations of pressure, saturation and temperature along the lapse of time, and the production amount of the shale gas in practical level can be achieved by the variation amount about the drilling well.

As described above, when a flow model for multi phase fluid is created, the flow model is activated by inputting initial conditions to each cells, and then inputting external conditions to each cells. The initial conditions to be inputted to the cells are initial pressure, saturation and temperature. And the external conditions such as depth of the drilling well, bottom pressure or a production amount or injection amount of the fluids can be set.

The research team of the present invention paid attention to a necessity of a method for grasping and predicting a multi phase flow of fluid to be as similar to the multi phase flow of fluids in real reservoir by particularly reflecting a temperature condition therein for the fluid such as a shale gas which moves along the continuous fractures defined in the reservoir unlike fluids such as oil or natural gas of the conventional resources that move along the pores in the reservoir.

Thus, in the present invention, a discrete fracture network model adapted to modeling fluids that are moving through fractures is used, and oil, water, gas and heat are modeled as fluids moving along the fractures in the reservoir in order to copy the multi phase flow. The significance of the present invention is with provision of a base enabling the prediction of various real aspects of the shale gas flow in the reservoir by such a method.

Above all, the present invention enables a more reliable modeling for the fluid flow in the shale gas area by adding the heat flow to an existing discrete fracture network flow model since there has not been enough research carried out for the part of heat and temperature in the case of multi phase flow in the discrete fracture network flow model despite the temperature condition acting as a very important factor in the multi phase flow.

Additional features of the present disclosure will become apparent to those skilled in the art upon consideration of illustrative embodiments exemplifying the best mode of carrying out the disclosure as presently perceived.

The present invention is disclosed in the accompanying figures with reference to an embodiment to provide an example of various features and concepts related to the invention, not to limit the scope of the invention. One skilled in the relevant art will recognize that numerous variations and modifications can be made to the configurations described above without departing from the scope of the present invention. The true scope of protection of the invention is defined by the appended claims. 

1. A method for analyzing multi phase and heat flow of fluids in a reservoir, the method comprising: (a) specifying a fracture surface in a reservoir as multiple three-dimensional coordinates and specifying a three-dimensional coordinate of an intersecting line between the fracture surfaces to grasp a connectivity between the fractures surfaces; (b) forming a fracture network by partitioning the fracture surface into multiple cells; (c) generating a flow model of each of at least two fluids among oil, water, and gas that move in the fracture network, and heat in a mathematical equation having main variables of pressure, saturation and temperature in the cell; and (d) grasping the flow of the fluids and heat within the fracture network in the reservoir layer by specifying initial values of the pressure, saturation and temperature in each of the cells, and specifying production condition or injection condition of a drilling well connected to the fracture network to operate the flow model.
 2. The method of claim 1, wherein, when the fracture surface is partitioned into the multiple cells, the intersecting line between the fracture surfaces becomes a boundary line between the cells.
 3. The method of claim 1, wherein the flow model of each of the fluids in the fracture network is a non-linear partial differential equation which is represented by a sum of flow amounts of the respective fluids moving along a longitudinal direction and a width direction of the cell within the cell.
 4. The method of claim 1, wherein the fluid comprises oil and water, wherein the flow model of oil is represented by Equation 1 below, and wherein the flow model of water is represented by Equation 2 below. $\begin{matrix} {{{{\frac{\partial}{\partial L}\left\lbrack {T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; W} + Q_{o} + Q_{Io}} = {\left( {\varphi \; V_{f}} \right)\frac{\partial}{\partial t}\left( \frac{S_{o}}{B_{o}} \right)}} & {{Equation}\mspace{14mu} 1} \\ {{{{{\frac{\partial}{\partial L}\left\lbrack {T_{wL}\frac{\partial}{\partial L}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {T_{wW}\frac{\partial}{\partial W}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; W} + Q_{w} + Q_{Iw}} = {\left( {\varphi \; V_{f}} \right)\frac{\partial}{\partial t}\left( \frac{S_{w}}{B_{w}} \right)}},} & {{Equation}\mspace{14mu} 2} \end{matrix}$ where L is a cell length, W is a cell width, T_(L) is a longitudinal conductivity of fluid in cell, T_(W) is width directional conductivity of fluid in cell, p is a fluid pressure in cell, γ is a cell density, z is a cell depth, Q is an amount of fluid produced in or injected from the drilling well, Q_(I) is a flow amount of fluid on intersecting line between fracture surfaces, φ is a cell porosity, V_(f) is a volume of fracture surface, t is time, S is a saturation rate of fluid in cell, B is a volume coefficient (volume) of fluid in cell, subscript o represents oil, and subscript w represents water.
 5. The method of claim 1, wherein the fluid comprises gas, and wherein the flow model of gas is represented by summing up the flow amounts of gas in state of gas and in state of gas being dissolved in oil.
 6. The method in claim 5, wherein the flow model of the gas is represented by Equation 3 below: $\begin{matrix} {{{{{\frac{\partial}{\partial L}\left\lbrack {{T_{gL}\frac{\partial}{\partial L}\left( {p_{g} - {\gamma_{g}z}} \right)} + {R_{so}T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)}} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {{T_{gW}\frac{\partial}{\partial W}\left( {p_{g} - {\gamma_{g}z}} \right)} + {R_{so}T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)}} \right\rbrack}\Delta \; W} + Q_{g} + {R_{so}Q_{o}} + Q_{Ig} + {R_{so}Q_{Io}}} = {\left( {\varphi \; V_{f}} \right)\frac{\partial}{\partial t}\left( {\frac{S_{g}}{B_{g}} + \frac{R_{so}S_{o}}{B_{o}}} \right)}},} & {{Equation}\mspace{14mu} 3} \end{matrix}$ where L is a cell length, W is a cell width, T_(L) is a longitudinal conductivity of fluid in cell, T_(W) is a width directional conductivity in cell, p is a pressure of fluid in cell, γ is a cell density, z is a cell depth, R_(s) is an amount of gas dissolved in another fluid, Q is an amount of fluid produced in or injected from the drilling well, Q_(I) is a flow amount on the intersecting line between the fractures surfaces, φ is a cell porosity, V_(f) is a volume of the fracture surface, t is time, S is a saturation rate of fluid in cell, B is a volume coefficient of fluid (volume) in cell, subscript o represents oil, and subscript g represents gas.
 7. The method of claim 1, wherein the flow model of heat in the fracture network is represented by a sum of temperature variation amounts in a longitudinal direction and a width direction in the cell, and a value obtained by multiplying enthalpy by the flow models of the fluids.
 8. The method in claim 7, wherein the flow model of heat is represented by Equation 4 below: $\begin{matrix} {{{{\frac{\partial}{\partial L}\left\lbrack {H_{o}T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {H_{w}T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)} \right\rbrack}\Delta \; W} + {{\frac{\partial}{\partial L}\left\lbrack {{H_{g}T_{gL}\frac{\partial}{\partial L}\left( {p_{g} - {\gamma_{g}z}} \right)} + {H_{g}T_{oL}\frac{\partial}{\partial L}\left( {p_{o} - {\gamma_{o}z}} \right)}} \right\rbrack}\Delta \; {L++}{\frac{\partial}{\partial W}\left\lbrack {{H_{g}T_{gW}\frac{\partial}{\partial W}\left( {p_{g} - {\gamma_{g}z}} \right)} + {H_{o}R_{so}T_{oW}\frac{\partial}{\partial W}\left( {p_{o} - {\gamma_{o}z}} \right)}} \right\rbrack}\Delta \; W} + {{\frac{\partial}{\partial L}\left\lbrack {H_{w}T_{wL}\frac{\partial}{\partial L}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; L} + {{\frac{\partial}{\partial W}\left\lbrack {H_{w}T_{wW}\frac{\partial}{\partial W}\left( {p_{w} - {\gamma_{w}z}} \right)} \right\rbrack}\Delta \; W} + {H_{o}\left( {Q_{o} + Q_{Io}} \right)} + {H_{g}\left\lbrack {Q_{g} + Q_{Ig} + {R_{so}\left( {Q_{o} + Q_{Io}} \right)}} \right\rbrack} + {H_{w}\left( {Q_{w} + Q_{Iw}} \right)} + {\frac{\partial}{\partial L}\left( {T_{h}\frac{\partial\theta}{\partial L}} \right)} + {\frac{\partial}{\partial W}\left( {T_{h}\frac{\partial\theta}{\partial W}} \right)}} = {\varphi \; V_{f}\frac{\partial}{\partial t}(U)}} & {{Equation}\mspace{14mu} 4} \end{matrix}$ where L is a cell length, W is a cell width, H is enthalpy of fluid in cell, T_(L) is conductivity of fluid in a longitudinal direction in cell, T_(W) is conductivity of fluid in a width direction in cell, p is a fluid pressure in cell, γ is a cell density, z is a cell depth, R_(s) is an amount of gas dissolved in another fluid, Q is an amount of fluid produced in or injected from the drilling well, Q_(I) is a flow amount of fluid on intersecting line between fracture surfaces, φ is a cell porosity, V_(f) is a volume of fracture surface, t is time, S is a saturation rate of fluid in cell, B is a volume coefficient (volume) in cell, θ is a temperature in cell, U is internal energy, subscript o represents oil, subscript w represents water, subscript g represents gas, and subscript h represents heat.
 9. A computer-readable recording media in which the method of any one of claims 1-8 is implemented by a program to be executed by a computer. 